The packages are automatically loaded using pacman. The
reported .html was last ran on the system: x86_64-apple-darwin17.0 and R
version: R version 4.2.1 (2022-06-23)
First, for the Stage 1 registered report some ICC
distributions are simulated, aiming for a bimodial distribution to pick
from and some specific categories (e.g., low, high, medium). Since
rnorm will simulate beyond what is preferred, subseting
< 1.00 and > -(.15).
# three diff ICC groups
# Elliott et al & Noble et al have largely reported a concentration around .3-.5. Here, simulate
# based on four levels, where min > -.15 and max is < 1.00
low_neg <- rnorm(10000, mean = 0, sd = 0.05)
low <- rnorm(10000, mean = 0.15, sd = 0.05) %>% subset(. < 1 & . > -.15)
mid <- rnorm(10000, mean = 0.45, sd = 0.1) %>% subset(. < 1 & . > -.15)
high <- rnorm(10000, mean = 0.7, sd = 0.2) %>% subset(. < 1 & . > -.15)
values <- c(low_neg, low, mid,high)plotting histogram for all values
hist(values, breaks = 10, col = "gray", xlab = "Values", main = "Bimodal Distribution of ICCs")Here, permutations are based on FWHM = five options, Motion = six options, Task Models = three options, Task Contrasts = four options, across three studies total 360 options that will be the “final” simulated dataset.
fwhm_ops = c(4,5,6,7,8) # five FWHM
motion_ops = c("None", "Regress_TransRot","Regress_Derivatives","Regress_1st8aCompCor",
"Censor_mfd","Excl_mFD") # fd cutoff, => .9
taskmod_ops = c("Cue_CueDur","Cue_CueFixDur","Fix_FixDur")
taskcontrast_ops = c("BigWin_v_Neut","BigWin_v_Base","SmWin_v_Neut","SmWin_v_Base")
studies = c("MLS","ABCD","AHRB")
# permutations of the variables & setting up for loop
var_permutations <- as.matrix(expand.grid(
study = studies, fwhm = fwhm_ops, motion = motion_ops,
modeltype = taskmod_ops, contrast = taskcontrast_ops))Creating empty DF to store info in
df = data.frame(study = character(), fwhm = numeric(),
motion = character(), contrast = character(), modeltype = character(),
med_icc = numeric())Here, using an abstract method (i.e., several if/else rules) to specify from which category of ICCs to sample from for specific modeling options. Note, this doesnt result in a covariance within studies, so the within-subj variance is likely to be random thus zero.
set.seed(1000)
# simulate across study
for (row in 1:nrow(var_permutations)) {
tmp_fwhm = var_permutations[row,"fwhm"]
tmp_motion = var_permutations[row,"motion"]
tmp_contrast = var_permutations[row,"contrast"]
tmp_model = var_permutations[row,"modeltype"]
tmp_study = var_permutations[row,"study"]
df[row,"fwhm"] <- tmp_fwhm
df[row,"motion"] <- tmp_motion
df[row,"contrast"] <- tmp_contrast
df[row,"modeltype"] <- tmp_model
df[row,"study"] <- tmp_study
if (tmp_fwhm %in% c(7,8)) {
if (tmp_motion %in% c("Regress_TransRot", "Regress_Derivatives","Regress_1st8aCompCor")) {
if (tmp_model %in% c("Cue_CueDur","Cue_CueFixDur")) {
if (tmp_contrast=="BigWin_v_Neut") {
df[row,"med_icc"] <- sample(x = high, size = 1,replace = FALSE)
} else {
df[row,"med_icc"] <- sample(x = mid, size = 1,replace = FALSE)
}
} else {
df[row,"med_icc"] <- sample(x = values, size = 1,replace = FALSE)
}
} else {
df[row,"med_icc"] <- sample(x = values, size = 1,replace = FALSE)
}
} else {
if (tmp_motion %in% c("Regress_TransRot", "Regress_Derivatives","Regress_1st8aCompCor")) {
if (tmp_model %in% c("Cue_CueDur","Cue_CueFixDur")) {
if (tmp_contrast=="BigWin_v_Neut") {
df[row,"med_icc"] <- sample(x = mid, size = 1,replace = FALSE)
} else {
df[row,"med_icc"] <- sample(x = low, size = 1,replace = FALSE)
}
} else {
df[row,"med_icc"] <- sample(x = low, size = 1,replace = FALSE)
}
} else {
df[row,"med_icc"] <- sample(x = low_neg, size = 1,replace = FALSE)
}
}
}releving and specifying factors so things are interpretable in the heiarchical models.
# level/relabel
df$fwhm_r <- as.numeric(df$fwhm)
df$motion = relevel(as.factor(df$motion), ref = "None",
levels = c("None","Censor_mfd","Excl_mFD",
"Regress_1st8aCompCor","Regress_Derivatives",
"Regress_TransRot"))
df$modeltype = relevel(as.factor(df$modeltype), ref = "Cue_CueDur",
levels = c("Cue_CueDur","Cue_CuFixDur","Fix_FixDur"))
df$contrast = relevel(as.factor(df$contrast), ref = "BigWin_v_Base",
levels = c("BigWin_v_Base","BigWin_v_Neut","SmWin_v_Base","SmWin_v_Neut"))Below, running the steps to summarize the different ICCs from the model combinations 360 for each study
Plotting overall and for each of [four] categories
icc_dist = df %>%
ggplot(aes(x = med_icc)) +
geom_histogram(bins = 50, fill = "white", colour = "black") +
ggtitle("Overall Distribution") +
theme_minimal()
fwhm_dist = df %>%
ggplot(aes(x = med_icc)) +
geom_histogram(bins = 50, fill = "white", colour = "black") +
facet_wrap(~fwhm) +
ggtitle("Distribution by FWHM category") +
theme_minimal()
motion_dist = df %>%
ggplot(aes(x = med_icc)) +
geom_histogram(bins = 50, fill = "white", colour = "black") +
facet_wrap(~motion) +
ggtitle("Distribution by Motion category") +
theme_minimal()
contrast_dist = df %>%
ggplot(aes(x = med_icc)) +
geom_histogram(bins = 50, fill = "white", colour = "black") +
facet_wrap(~contrast) +
ggtitle("Distribution by Contrast category") +
theme_minimal()
modeltype_dist = df %>%
ggplot(aes(x = med_icc)) +
geom_histogram(bins = 50, fill = "white", colour = "black") +
facet_wrap(~modeltype) +
ggtitle("Distribution by Model Type category") +
theme_minimal()
icc_dist fwhm_distmotion_distcontrast_distmodeltype_dist
Creating rainclouds via ggrain
fwhm_rg = df %>% ggplot(aes(x = fwhm, y = med_icc, fill = fwhm, color = fwhm)) +
geom_rain(alpha = .5, rain.side = 'l',
boxplot.args = list(color = "black", outlier.shape = NA),
boxplot.args.pos = list(
position = ggpp::position_dodgenudge(x = .1), width = 0.1
)) +
theme_classic() +
ggtitle("Distribution by FWHN category") +
scale_fill_brewer(palette = 'Dark2') +
scale_color_brewer(palette = 'Dark2') +
guides(fill = 'none', color = 'none') +
theme(text = element_text(family = "Times New Roman"),
axis.text = element_text(size = 10),
axis.title = element_text(size = 10),
legend.text = element_text(size = 10),
legend.title = element_text(size = 10),
plot.title = element_text(size = 16))
motion_rg = df %>% ggplot(aes(x = motion, y = med_icc, fill = motion, color = motion)) +
geom_rain(alpha = .5, rain.side = 'l',
boxplot.args = list(color = "black", outlier.shape = NA),
boxplot.args.pos = list(
position = ggpp::position_dodgenudge(x = .1), width = 0.1
)) +
theme_classic() +
ggtitle("Distribution by Motion category") +
scale_fill_brewer(palette = 'Dark2') +
scale_color_brewer(palette = 'Dark2') +
guides(fill = 'none', color = 'none') +
theme(text = element_text(family = "Times New Roman"),
axis.text = element_text(size = 10),
axis.title = element_text(size = 10),
legend.text = element_text(size = 10),
legend.title = element_text(size = 10),
plot.title = element_text(size = 16))
modeltype_rg = df %>% ggplot(aes(x = modeltype, y = med_icc, fill = modeltype, color = modeltype)) +
geom_rain(alpha = .5, rain.side = 'l',
boxplot.args = list(color = "black", outlier.shape = NA),
boxplot.args.pos = list(
position = ggpp::position_dodgenudge(x = .1), width = 0.1
)) +
theme_classic() +
ggtitle("Distribution by Model Type category") +
scale_fill_brewer(palette = 'Dark2') +
scale_color_brewer(palette = 'Dark2') +
guides(fill = 'none', color = 'none') +
theme(text = element_text(family = "Times New Roman"),
axis.text = element_text(size = 10),
axis.title = element_text(size = 10),
legend.text = element_text(size = 10),
legend.title = element_text(size = 10),
plot.title = element_text(size = 16))
contrast_rg = df %>% ggplot(aes(x = contrast, y = med_icc, fill = contrast, color = contrast)) +
geom_rain(alpha = .5, rain.side = 'l',
boxplot.args = list(color = "black", outlier.shape = NA),
boxplot.args.pos = list(
position = ggpp::position_dodgenudge(x = .1), width = 0.1
)) +
theme_classic() +
ggtitle("Distribution by Contrast category") +
scale_fill_brewer(palette = 'Dark2') +
scale_color_brewer(palette = 'Dark2') +
guides(fill = 'none', color = 'none') +
theme(text = element_text(family = "Times New Roman"),
axis.text = element_text(size = 10),
axis.title = element_text(size = 10),
legend.text = element_text(size = 10),
legend.title = element_text(size = 10),
plot.title = element_text(size = 16))
fwhm_rgmotion_rgmodeltype_rgcontrast_rgsummary_icc = df %>%
group_by(fwhm, motion, modeltype, contrast) %>%
summarise(across(("med_icc"), list(min = min, max = max, mean = mean, median = median))) %>%
as.data.frame()## `summarise()` has grouped output by 'fwhm', 'motion', 'modeltype'. You can
## override using the `.groups` argument.
summary_icc %>%
arrange(desc(med_icc_min)) %>%
slice(1:10)## fwhm motion modeltype contrast med_icc_min
## 1 7 Regress_1st8aCompCor Cue_CueFixDur BigWin_v_Neut 0.7562752
## 2 7 Regress_TransRot Cue_CueDur BigWin_v_Neut 0.7342187
## 3 7 Regress_1st8aCompCor Cue_CueDur BigWin_v_Neut 0.6781959
## 4 8 Regress_1st8aCompCor Cue_CueDur BigWin_v_Neut 0.6599825
## 5 8 Regress_TransRot Cue_CueFixDur BigWin_v_Neut 0.6202459
## 6 8 Regress_TransRot Cue_CueDur BigWin_v_Neut 0.6006907
## 7 7 Regress_Derivatives Cue_CueDur BigWin_v_Neut 0.5129165
## 8 7 Regress_Derivatives Cue_CueDur SmWin_v_Base 0.5122227
## 9 8 Regress_TransRot Cue_CueDur SmWin_v_Base 0.4910063
## 10 7 Censor_mfd Cue_CueFixDur BigWin_v_Base 0.4752819
## med_icc_max med_icc_mean med_icc_median
## 1 0.9152130 0.8166406 0.7784335
## 2 0.8155691 0.7797988 0.7896087
## 3 0.8508815 0.7753146 0.7968663
## 4 0.8927606 0.7572947 0.7191411
## 5 0.8024149 0.7262381 0.7560534
## 6 0.7947373 0.6947145 0.6887155
## 7 0.6764460 0.5986463 0.6065764
## 8 0.6208878 0.5716596 0.5818682
## 9 0.5377776 0.5122792 0.5080539
## 10 0.7166981 0.6078405 0.6315416
summary_icc %>%
arrange(desc(med_icc_max)) %>%
slice(1:10)## fwhm motion modeltype contrast med_icc_min
## 1 7 Censor_mfd Fix_FixDur SmWin_v_Base 0.02215935
## 2 7 Censor_mfd Fix_FixDur BigWin_v_Base 0.02659101
## 3 8 Excl_mFD Cue_CueFixDur SmWin_v_Neut 0.10338313
## 4 7 Excl_mFD Cue_CueDur BigWin_v_Base 0.08332403
## 5 7 Regress_TransRot Fix_FixDur SmWin_v_Neut 0.41721161
## 6 7 Regress_1st8aCompCor Cue_CueFixDur BigWin_v_Neut 0.75627515
## 7 8 Excl_mFD Cue_CueDur BigWin_v_Neut 0.12289872
## 8 8 None Cue_CueDur SmWin_v_Neut -0.04767462
## 9 8 Regress_1st8aCompCor Cue_CueDur BigWin_v_Neut 0.65998255
## 10 7 None Fix_FixDur SmWin_v_Neut -0.05309287
## med_icc_max med_icc_mean med_icc_median
## 1 0.9947495 0.6162410 0.8318142
## 2 0.9640912 0.4915768 0.4840481
## 3 0.9435615 0.3880643 0.1172483
## 4 0.9398068 0.4413431 0.3008984
## 5 0.9364607 0.6109669 0.4792284
## 6 0.9152130 0.8166406 0.7784335
## 7 0.9014496 0.5025672 0.4833531
## 8 0.8950570 0.4689903 0.5595883
## 9 0.8927606 0.7572947 0.7191411
## 10 0.8536519 0.3142251 0.1421162
summary_icc %>%
arrange(desc(med_icc_mean)) %>%
slice(1:10)## fwhm motion modeltype contrast med_icc_min
## 1 7 Regress_1st8aCompCor Cue_CueFixDur BigWin_v_Neut 0.75627515
## 2 7 Regress_TransRot Cue_CueDur BigWin_v_Neut 0.73421866
## 3 7 Regress_1st8aCompCor Cue_CueDur BigWin_v_Neut 0.67819590
## 4 8 Regress_1st8aCompCor Cue_CueDur BigWin_v_Neut 0.65998255
## 5 8 Regress_TransRot Cue_CueFixDur BigWin_v_Neut 0.62024593
## 6 8 Regress_TransRot Cue_CueDur BigWin_v_Neut 0.60069068
## 7 8 Regress_1st8aCompCor Cue_CueFixDur BigWin_v_Neut 0.42787099
## 8 8 Regress_Derivatives Cue_CueFixDur BigWin_v_Neut 0.44973463
## 9 7 Regress_Derivatives Cue_CueFixDur BigWin_v_Neut 0.44215137
## 10 7 Censor_mfd Fix_FixDur SmWin_v_Base 0.02215935
## med_icc_max med_icc_mean med_icc_median
## 1 0.9152130 0.8166406 0.7784335
## 2 0.8155691 0.7797988 0.7896087
## 3 0.8508815 0.7753146 0.7968663
## 4 0.8927606 0.7572947 0.7191411
## 5 0.8024149 0.7262381 0.7560534
## 6 0.7947373 0.6947145 0.6887155
## 7 0.8232779 0.6800934 0.7891313
## 8 0.8360832 0.6592579 0.6919558
## 9 0.7777469 0.6316792 0.6751393
## 10 0.9947495 0.6162410 0.8318142
Below using lmer to fit HLM model. One model has random
intercept (study), second model has random intercept (study) and random
slope (FWHM). Reporting using tab_model
Level 1 model
\(MedianICC_{ij} = \beta_{0_j} + \beta_{1_j} \times \text{fwhm}_{ij} + \beta_{2_j} \times \text{motion}_{ij} + \beta_{3_j} \times \text{modeltype}_{ij} + \beta_{4_j} \times \text{contrast}_{ij} + \epsilon_{ij}\)
Level 2 model
\(\beta_{0_j} = \gamma_{00} + u_{0_j}\)
\(\beta_{1_j} = \gamma_{10} + u_{1_j}\)
\(\beta_{2_j} = \gamma_{20} + u_{2_j}\)
\(\beta_{3_j} = \gamma_{30} + u_{3_j}\)
\(\beta_{4_j} = \gamma_{40} + u_{4_j}\)
`Notes:
med_icc_ij is the outcome for the ith model in the jth sample β0_j is the intercept for the jth sample fwhm_r_ij, motion_ij, modeltype_ij, and contrast_ij are the predictor variables for the ith model in the jth sample β1_j, β2_j, β3_j, and β4_j: corresponding coefficients for the jth sample ε_ij is the within-sample error term
γ00, γ10, γ20, γ30, and γ40 are the fixed effects intercept and slopes for the predictors; u0_j, u1_j, u2_j, u3_j, and u4_j are the sample-level residuals (random effects) for the intercept and the four slopes.`
model = lmer(med_icc ~ fwhm_r + motion + modeltype + contrast + (1 | study), data = df)
model_rndslp = lmer(med_icc ~ fwhm_r + motion + modeltype + contrast + (1+fwhm_r| study), data = df)## boundary (singular) fit: see help('isSingular')
tab_model(model, model_rndslp)| med_icc | med_icc | |||||
|---|---|---|---|---|---|---|
| Predictors | Estimates | CI | p | Estimates | CI | p |
| (Intercept) | -0.35 | -0.41 – -0.28 | <0.001 | -0.35 | -0.42 – -0.27 | <0.001 |
| fwhm r | 0.08 | 0.07 – 0.09 | <0.001 | 0.08 | 0.07 – 0.09 | <0.001 |
| motion [Censor_mfd] | -0.02 | -0.06 – 0.02 | 0.382 | -0.02 | -0.06 – 0.02 | 0.381 |
| motion [Excl_mFD] | 0.01 | -0.03 – 0.05 | 0.667 | 0.01 | -0.03 – 0.05 | 0.667 |
|
motion [Regress_1st8aCompCor] |
0.16 | 0.12 – 0.20 | <0.001 | 0.16 | 0.12 – 0.20 | <0.001 |
|
motion [Regress_Derivatives] |
0.15 | 0.11 – 0.18 | <0.001 | 0.15 | 0.11 – 0.18 | <0.001 |
| motion [Regress_TransRot] | 0.17 | 0.13 – 0.21 | <0.001 | 0.17 | 0.13 – 0.21 | <0.001 |
| modeltype [Cue_CueFixDur] | -0.01 | -0.04 – 0.02 | 0.432 | -0.01 | -0.04 – 0.02 | 0.431 |
| modeltype [Fix_FixDur] | -0.08 | -0.10 – -0.05 | <0.001 | -0.08 | -0.10 – -0.05 | <0.001 |
| contrast [BigWin_v_Neut] | 0.08 | 0.05 – 0.12 | <0.001 | 0.08 | 0.05 – 0.12 | <0.001 |
| contrast [SmWin_v_Base] | 0.00 | -0.03 – 0.03 | 0.991 | 0.00 | -0.03 – 0.03 | 0.991 |
| contrast [SmWin_v_Neut] | -0.01 | -0.04 – 0.03 | 0.721 | -0.01 | -0.04 – 0.03 | 0.720 |
| Random Effects | ||||||
| σ2 | 0.04 | 0.04 | ||||
| τ00 | 0.00 study | 0.00 study | ||||
| τ11 | 0.00 study.fwhm_r | |||||
| ρ01 | -1.00 study | |||||
| ICC | 0.00 | |||||
| N | 3 study | 3 study | ||||
| Observations | 1080 | 1080 | ||||
| Marginal R2 / Conditional R2 | 0.383 / 0.383 | 0.384 / NA | ||||
Using emmeans to control type I error rate of controls via Tukey’s Honest Significant Test
# Model comparisons
mod1_lmer_obs = as.numeric(summary(model)[3]$devcomp$dims[1])
emm_options(lmerTest.limit = mod1_lmer_obs)
# contrast contrast
mod1_thsd = emmeans(object = model, specs = pairwise ~ contrast, pbkrtest.limit = mod1_lmer_obs)
mod1_thsd$contrasts %>% summary(infer = TRUE) %>%
kbl(digits = 4, caption = "Contrast Tukey HSD: Contrast") %>% kable_styling()| contrast | estimate | SE | df | lower.CL | upper.CL | t.ratio | p.value |
|---|---|---|---|---|---|---|---|
| BigWin_v_Base - BigWin_v_Neut | -0.0846 | 0.0163 | 1066 | -0.1265 | -0.0426 | -5.1852 | 0.0000 |
| BigWin_v_Base - SmWin_v_Base | -0.0002 | 0.0163 | 1066 | -0.0421 | 0.0418 | -0.0110 | 1.0000 |
| BigWin_v_Base - SmWin_v_Neut | 0.0058 | 0.0163 | 1066 | -0.0361 | 0.0478 | 0.3578 | 0.9843 |
| BigWin_v_Neut - SmWin_v_Base | 0.0844 | 0.0163 | 1066 | 0.0424 | 0.1264 | 5.1742 | 0.0000 |
| BigWin_v_Neut - SmWin_v_Neut | 0.0904 | 0.0163 | 1066 | 0.0484 | 0.1324 | 5.5430 | 0.0000 |
| SmWin_v_Base - SmWin_v_Neut | 0.0060 | 0.0163 | 1066 | -0.0360 | 0.0480 | 0.3688 | 0.9829 |
# contrast motion
mod1_thsd = emmeans(object = model, specs = pairwise ~ motion, pbkrtest.limit = mod1_lmer_obs)
mod1_thsd$contrasts %>% summary(infer = TRUE) %>%
kbl(digits = 4, caption = "Contrast Tukey HSD: Motion Corr") %>% kable_styling()| contrast | estimate | SE | df | lower.CL | upper.CL | t.ratio | p.value |
|---|---|---|---|---|---|---|---|
| None - Censor_mfd | 0.0175 | 0.02 | 1066 | -0.0396 | 0.0745 | 0.8738 | 0.9527 |
| None - Excl_mFD | -0.0086 | 0.02 | 1066 | -0.0656 | 0.0484 | -0.4297 | 0.9981 |
| None - Regress_1st8aCompCor | -0.1605 | 0.02 | 1066 | -0.2175 | -0.1035 | -8.0352 | 0.0000 |
| None - Regress_Derivatives | -0.1455 | 0.02 | 1066 | -0.2026 | -0.0885 | -7.2862 | 0.0000 |
| None - Regress_TransRot | -0.1666 | 0.02 | 1066 | -0.2236 | -0.1095 | -8.3391 | 0.0000 |
| Censor_mfd - Excl_mFD | -0.0260 | 0.02 | 1066 | -0.0831 | 0.0310 | -1.3035 | 0.7832 |
| Censor_mfd - Regress_1st8aCompCor | -0.1780 | 0.02 | 1066 | -0.2350 | -0.1209 | -8.9089 | 0.0000 |
| Censor_mfd - Regress_Derivatives | -0.1630 | 0.02 | 1066 | -0.2200 | -0.1060 | -8.1599 | 0.0000 |
| Censor_mfd - Regress_TransRot | -0.1840 | 0.02 | 1066 | -0.2410 | -0.1270 | -9.2128 | 0.0000 |
| Excl_mFD - Regress_1st8aCompCor | -0.1519 | 0.02 | 1066 | -0.2089 | -0.0949 | -7.6054 | 0.0000 |
| Excl_mFD - Regress_Derivatives | -0.1370 | 0.02 | 1066 | -0.1940 | -0.0799 | -6.8564 | 0.0000 |
| Excl_mFD - Regress_TransRot | -0.1580 | 0.02 | 1066 | -0.2150 | -0.1010 | -7.9093 | 0.0000 |
| Regress_1st8aCompCor - Regress_Derivatives | 0.0150 | 0.02 | 1066 | -0.0421 | 0.0720 | 0.7490 | 0.9756 |
| Regress_1st8aCompCor - Regress_TransRot | -0.0061 | 0.02 | 1066 | -0.0631 | 0.0510 | -0.3039 | 0.9997 |
| Regress_Derivatives - Regress_TransRot | -0.0210 | 0.02 | 1066 | -0.0781 | 0.0360 | -1.0529 | 0.8997 |
# contrast modeltype
mod1_thsd = emmeans(object = model, specs = pairwise ~ modeltype, pbkrtest.limit = mod1_lmer_obs)
mod1_thsd$contrasts %>% summary(infer = TRUE) %>%
kbl(digits = 4, caption = "Contrast Tukey HSD: Model Type ") %>% kable_styling()| contrast | estimate | SE | df | lower.CL | upper.CL | t.ratio | p.value |
|---|---|---|---|---|---|---|---|
| Cue_CueDur - Cue_CueFixDur | 0.0111 | 0.0141 | 1066 | -0.0220 | 0.0442 | 0.7858 | 0.7118 |
| Cue_CueDur - Fix_FixDur | 0.0764 | 0.0141 | 1066 | 0.0433 | 0.1096 | 5.4101 | 0.0000 |
| Cue_CueFixDur - Fix_FixDur | 0.0653 | 0.0141 | 1066 | 0.0322 | 0.0985 | 4.6242 | 0.0000 |
Creating data in a format that is compatible with specr. Needs: estimate (i.e., ICC), std.error, conf.high, conf.low.
# create specr plot for med_icc averages for
# first, combine independent model vars into string to create average for each model type
df$model_type <- paste(df$fwhm,df$motion,df$contrast,df$modeltype,sep = "-")
# calculate the avg estimate of ICC across study, standard error and +/- 95% confidence interval. In complete version
df_summ <- df %>%
group_by(model_type) %>%
summarise(estimate = mean(med_icc), std.error = sd(med_icc)/sqrt(length(med_icc))) %>%
mutate(conf.low = estimate - 1.96 * std.error, conf.high = estimate + 1.96 * std.error) %>%
separate(col = model_type, into = c("fwhm","motion","contrast","modeltype"), sep = "-", remove = FALSE)creating combined panel 1 and panel 2 for all model permutations first.
plot_a = plot_curve(df = df_summ, ci = TRUE, desc = FALSE, legend = FALSE, null = 0)
plot_b <- plot_choices(df = df_summ, choices = c("fwhm", "motion","contrast","modeltype"), desc = F, null = 0) +
labs(y = "Variables", x = "Ordered Specification Curve \n ICC coefficient")
cowplot::plot_grid(plot_a, plot_b, ncol = 1, align = "v", axis = 'tblr',
labels = c('A', 'B'), rel_heights = c(1, 2),
label_fontfamily = "Times", label_size = 12)Creating model that is subset to visualize reliability for the top (>75th) and bottom (< 25th) quartile
# get 75th/25th qunatiles
top_75q = as.numeric(quantile(df_summ$estimate, .75))
bot_25q = as.numeric(quantile(df_summ$estimate, .25))
df_summ_subset = df_summ %>%
filter(estimate < bot_25q | estimate > top_75q)
# subset plots
plot_a_sub = plot_curve(df = df_summ_subset, ci = TRUE, desc = FALSE, legend = FALSE, null = 0)
plot_b_sub <- plot_choices(df = df_summ_subset, choices = c("fwhm", "motion","contrast","modeltype"), desc = F, null = 0) +
labs(y = "Variables", x = "Ordered Specification Curve \n ICC coefficient")
cowplot::plot_grid(plot_a_sub, plot_b_sub, ncol = 1, align = "v", axis = 'tblr',
labels = c('A', 'B'), rel_heights = c(1, 2),
label_fontfamily = "Times", label_size = 12)The ability to estimate the BOLD signal in the timeseries, efficiency plays a significant role. See issues discussed by series of videos from Jeanette Mumford. The MID task, has historically different ways of defining the “Anticipation” phase. Previously, it started by modeling the fixation phase/duration and later some have used the cue/duration. However, the issues have been rarely cleared up in the literature as it is challenging to decipher from the information provided in papers what is defined as the anticipation phase and duration. Hence, up to three models can be created in the design matrix.
Below, using the efficiency calculation tutorial from Dr Mumford to estimate the effieiency of the three different models that can be estimated in this version of the MID task.
First, the behavior data is pulled from the events files for each subject and each run. Some subjects are excluded because they lack hit/miss trials which causes an error in the step using neuRosim.
# location for subjects run behavior data
# update path as needed
dir = "/Users/michaeldemidenko/Desktop/Academia/UM/2_AHRB/Projects/BIDS/Beh/MID_W1/"
# getting list of file names
beh_files = list.files(dir)
# removing subjects that dont have some Hit/Miss for outcome, as there are zero values for complete calculation
beh_files_ct = beh_files[-c(12,19,34,39,46,120,140)]
# sub-07_ses-1_task-mid_run-02_events.tsv
# sub-100_ses-1_task-mid_run-01_events.tsv
# sub-107_ses-1_task-mid_run-02_events.tsv
# sub-12_ses-1_task-mid_run-01_events.tsv
# sub-15_ses-1_task-mid_run-02_events.tsv
# sub-52_ses-1_task-mid_run-02_events.tsv
# sub-62_ses-1_task-mid_run-02_events.tsvUsing for loop to create the the data based on the behavioral onsets/durations for efficiency results. Three different models are fit:
design_CueMod: Cue Onset + Cue Durationdesign_AntMod: Cue Onset + (Cue Duration + Fixation
Duration)design_FixMod: Fixation Onset + Fixation DurationWhile the focus is on the anticipatory phase, e.g. modeling 5 cue times + associated duration, the feedback phase is modeled, too. In other words, the total model includes 15 items.
# Null df to add to
df = NULL
# creating forloop using Jeanette Mumfords recommended code from Neurohack - https://github.com/jmumford/nhwEfficiency
# tailored to work with behavior data collected for MID task data
#set TR, vols, effect size (amplitude) and run and mri length
tr = .8
vols = 407
effect_size = 1
run_len = vols*tr
mri_len = seq(1,407,1)
for (f in beh_files_ct) {
# Combining directory and file path to open .csv
beh_dat = read.csv(paste(dir,f,sep=""),sep = '\t')
# subj/run labels to use in later df creation
str_vals =data.frame(str_split(f,'_'))
sub = str_vals[1,1]
run = str_vals[4,1]
# convert one of the outcome items foor Neutral condition so it is differently labeled for Hit/Miss
beh_dat <- beh_dat %>%
mutate(TRIAL_RESULT = if_else(TRIAL_RESULT!="No money at stake!",TRIAL_RESULT,
if_else(PROBE_HIT==1,TRIAL_RESULT,"No Money At Stake, MISS")))
#
### Anticipation
# Onsets
BW_ons_1 <- beh_dat$CUE_ONSET[beh_dat$TRIAL_TYPE=="LargeGain"]
SW_ons_1 <- beh_dat$CUE_ONSET[beh_dat$TRIAL_TYPE=="SmallGain"]
Neut_ons_1 <- beh_dat$CUE_ONSET[beh_dat$TRIAL_TYPE=="NoMoneyStake"]
BL_ons_1 <- beh_dat$CUE_ONSET[beh_dat$TRIAL_TYPE=="LargeLoss"]
SL_ons_1 <- beh_dat$CUE_ONSET[beh_dat$TRIAL_TYPE=="SmallLoss"]
BW_ons_2 <- beh_dat$FIXATION_ONSET[beh_dat$TRIAL_TYPE=="LargeGain"]
SW_ons_2 <- beh_dat$FIXATION_ONSET[beh_dat$TRIAL_TYPE=="SmallGain"]
Neut_ons_2 <- beh_dat$FIXATION_ONSET[beh_dat$TRIAL_TYPE=="NoMoneyStake"]
BL_ons_2 <- beh_dat$FIXATION_ONSET[beh_dat$TRIAL_TYPE=="LargeLoss"]
SL_ons_2 <- beh_dat$FIXATION_ONSET[beh_dat$TRIAL_TYPE=="SmallLoss"]
# durations
BW_dur1 <- beh_dat$CUE_DURATION[beh_dat$TRIAL_TYPE=="LargeGain"]
SW_dur1 <- beh_dat$CUE_DURATION[beh_dat$TRIAL_TYPE=="SmallGain"]
Neut_dur1 <- beh_dat$CUE_DURATION[beh_dat$TRIAL_TYPE=="NoMoneyStake"]
BL_dur1 <- beh_dat$CUE_DURATION[beh_dat$TRIAL_TYPE=="LargeLoss"]
SL_dur1 <- beh_dat$CUE_DURATION[beh_dat$TRIAL_TYPE=="SmallLoss"]
BW_dur2 <- beh_dat$FIXATION_DURATION[beh_dat$TRIAL_TYPE=="LargeGain"]
SW_dur2 <- beh_dat$FIXATION_DURATION[beh_dat$TRIAL_TYPE=="SmallGain"]
Neut_dur2 <- beh_dat$FIXATION_DURATION[beh_dat$TRIAL_TYPE=="NoMoneyStake"]
BL_dur2 <- beh_dat$FIXATION_DURATION[beh_dat$TRIAL_TYPE=="LargeLoss"]
SL_dur2 <- beh_dat$FIXATION_DURATION[beh_dat$TRIAL_TYPE=="SmallLoss"]
BW_dur12 <- beh_dat$CUE_DURATION[beh_dat$TRIAL_TYPE=="LargeGain"] + beh_dat$FIXATION_DURATION[beh_dat$TRIAL_TYPE=="LargeGain"]
SW_dur12 <- beh_dat$CUE_DURATION[beh_dat$TRIAL_TYPE=="SmallGain"] + beh_dat$FIXATION_DURATION[beh_dat$TRIAL_TYPE=="SmallGain"]
Neut_dur12 <- beh_dat$CUE_DURATION[beh_dat$TRIAL_TYPE=="NoMoneyStake"]+ beh_dat$FIXATION_DURATION[beh_dat$TRIAL_TYPE=="NoMoneyStake"]
BL_dur12 <- beh_dat$CUE_DURATION[beh_dat$TRIAL_TYPE=="LargeLoss"] + beh_dat$FIXATION_DURATION[beh_dat$TRIAL_TYPE=="LargeLoss"]
SL_dur12 <- beh_dat$CUE_DURATION[beh_dat$TRIAL_TYPE=="SmallLoss"] + beh_dat$FIXATION_DURATION[beh_dat$TRIAL_TYPE=="SmallLoss"]
### Feedback
# Onsets
BW_Hit_ons <- beh_dat$FEEDBACK_ONSET[beh_dat$TRIAL_RESULT=="You earn $5!"]
SW_Hit_ons <- beh_dat$FEEDBACK_ONSET[beh_dat$TRIAL_RESULT=="You earn $0.20!"]
Neut_Hit_ons <- beh_dat$FEEDBACK_ONSET[beh_dat$TRIAL_RESULT=="No money at stake!"]
BL_Hit_ons <- beh_dat$FEEDBACK_ONSET[beh_dat$TRIAL_RESULT=="You keep $5!"]
SL_Hit_ons <- beh_dat$FEEDBACK_ONSET[beh_dat$TRIAL_RESULT=="You keep $0.20!"]
BW_Miss_ons <- beh_dat$FEEDBACK_ONSET[beh_dat$TRIAL_RESULT=="You did not earn $5!"]
SW_Miss_ons <- beh_dat$FEEDBACK_ONSET[beh_dat$TRIAL_RESULT=="You did not earn $0.20!"]
Neut_Miss_ons <- beh_dat$FEEDBACK_ONSET[beh_dat$TRIAL_RESULT=="No Money At Stake, MISS"]
BL_Miss_ons <- beh_dat$FEEDBACK_ONSET[beh_dat$TRIAL_RESULT=="You lose $5!"]
SL_Miss_ons <- beh_dat$FEEDBACK_ONSET[beh_dat$TRIAL_RESULT=="You lose $0.20!"]
# durations
BW_Hit_dur <- beh_dat$FEEDBACK_DURATION[beh_dat$TRIAL_RESULT=="You earn $5!"]
SW_Hit_dur <- beh_dat$FEEDBACK_DURATION[beh_dat$TRIAL_RESULT=="You earn $0.20!"]
Neut_Hit_dur <- beh_dat$FEEDBACK_DURATION[beh_dat$TRIAL_RESULT=="No money at stake!"]
BL_Hit_dur <- beh_dat$FEEDBACK_DURATION[beh_dat$TRIAL_RESULT=="You keep $5!"]
SL_Hit_dur <- beh_dat$FEEDBACK_DURATION[beh_dat$TRIAL_RESULT=="You keep $0.20!"]
BW_Miss_dur <- beh_dat$FEEDBACK_DURATION[beh_dat$TRIAL_RESULT=="You did not earn $5!"]
SW_Miss_dur <- beh_dat$FEEDBACK_DURATION[beh_dat$TRIAL_RESULT=="You did not earn $0.20!"]
Neut_Miss_dur <- beh_dat$FEEDBACK_DURATION[beh_dat$TRIAL_RESULT=="No Money At Stake, MISS"]
BL_Miss_dur <- beh_dat$FEEDBACK_DURATION[beh_dat$TRIAL_RESULT=="You lose $5!"]
SL_Miss_dur <- beh_dat$FEEDBACK_DURATION[beh_dat$TRIAL_RESULT=="You lose $0.20!"]
## model 1 (cue only), creating a list of values for Onset and Duration
CueMod_ons = list(BW_ons_1, SW_ons_1,Neut_ons_1,BL_ons_1,SL_ons_1,
BW_Hit_ons, SW_Hit_ons, Neut_Hit_ons, BL_Hit_ons,SL_Hit_ons,
BW_Miss_ons,SW_Miss_ons, Neut_Miss_ons, BL_Miss_ons,SL_Miss_ons
)
CueMod_dur = list(BW_dur1,SW_dur1,Neut_dur1,BL_dur1,SL_dur1,
BW_Hit_dur, SW_Hit_dur, Neut_Hit_dur, BL_Hit_dur,SL_Hit_dur,
BW_Miss_dur,SW_Miss_dur, Neut_Miss_dur, BL_Miss_dur,SL_Miss_dur
)
# Using specifydesign to create the HRF Cue Model based on inputs:
# onsets of the specified model
# durations of the specificated model
# the list of effect size (or amplitudes) used for each stimulus in model
# Tr and total time to create the length of HRF using the convolved double-gamma
# condition names is a list of labels for the Input onset labels
Cue_Mod <- specifydesign(onsets = CueMod_ons,
durations = CueMod_dur,
effectsize =list(1,1,1,1,1,
1,1,1,1,1,
1,1,1,1,1),
TR = tr, totaltime = run_len, conv = "double-gamma",
cond.names = list("BW","SW","Neut","BL","SL",
"BWhit","SWhit","NEUThit","BLhit","SLhit",
"BWmss","SWmss","NEUTmss","BLmss","SLmss")
)
# Model 2 (fixation only)
## model 1
FixMod_ons = list(BW_ons_2, SW_ons_2,Neut_ons_2,BL_ons_2,SL_ons_2,
BW_Hit_ons, SW_Hit_ons, Neut_Hit_ons, BL_Hit_ons,SL_Hit_ons,
BW_Miss_ons,SW_Miss_ons, Neut_Miss_ons, BL_Miss_ons,SL_Miss_ons
)
FixMod_dur = list(BW_dur2,SW_dur2,Neut_dur2,BL_dur2,SL_dur2,
BW_Hit_dur, SW_Hit_dur, Neut_Hit_dur, BL_Hit_dur,SL_Hit_dur,
BW_Miss_dur,SW_Miss_dur, Neut_Miss_dur, BL_Miss_dur,SL_Miss_dur
)
Fix_Mod <- specifydesign(onsets = FixMod_ons,
durations = FixMod_dur,
effectsize =list(1,1,1,1,1,
1,1,1,1,1,
1,1,1,1,1),
TR = tr, totaltime = run_len, conv = "double-gamma",
cond.names = list("BW","SW","Neut","BL","SL",
"BWhit","SWhit","NEUThit","BLhit","SLhit",
"BWmss","SWmss","NEUTmss","BLmss","SLmss")
)
#Anticipation model cue onset & cue + fix dur
## model 1
AntMod_ons = list(BW_ons_1, SW_ons_1,Neut_ons_1,BL_ons_1,SL_ons_1,
BW_Hit_ons, SW_Hit_ons, Neut_Hit_ons, BL_Hit_ons,SL_Hit_ons,
BW_Miss_ons,SW_Miss_ons, Neut_Miss_ons, BL_Miss_ons,SL_Miss_ons
)
AntMod_dur = list(BW_dur12,SW_dur12,Neut_dur12,BL_dur12,SL_dur12,
BW_Hit_dur, SW_Hit_dur, Neut_Hit_dur, BL_Hit_dur,SL_Hit_dur,
BW_Miss_dur,SW_Miss_dur, Neut_Miss_dur, BL_Miss_dur,SL_Miss_dur
)
Ant_Mod <- specifydesign(onsets = AntMod_ons,
durations = AntMod_dur,
effectsize =list(1,1,1,1,1,
1,1,1,1,1,
1,1,1,1,1),
TR = tr, totaltime = run_len, conv = "double-gamma",
cond.names = list("BW","SW","Neut","BL","SL",
"BWhit","SWhit","NEUThit","BLhit","SLhit",
"BWmss","SWmss","NEUTmss","BLmss","SLmss")
)
design_CueMod <- data.frame("TR" = mri_len, Cue_Mod)
# create design matrix + test contrast
CueMod_mat = cbind("Inter" = rep(1,length(design_CueMod$BW)),
design_CueMod[2:16])
# Plot fixation model and create design matrix
design_FixMod <- data.frame("TR" = mri_len, Fix_Mod)
# create design matrix
FixMod_mat = cbind("Inter" = rep(1,length(design_FixMod$BW)),
design_FixMod[2:16])
# Plot anticipation model and create design matrix
design_AntMod <- data.frame("TR" = mri_len, Ant_Mod)
# create design matrix
AntMod_mat = cbind("Inter" = rep(1,length(design_AntMod$BW)),
design_AntMod[2:16])
# setting up contrasts to estimate efficiency
LGain = c(0,
1,0,-1,0,0, # large gain versus neutral
0,0,0,0,0,0,0,0,0,0
)
LGain_CueMod = 1/(t(LGain)%*%solve(t(CueMod_mat)%*%as.matrix(CueMod_mat))%*%LGain)
LGain_AntMod = 1/(t(LGain)%*%solve(t(AntMod_mat)%*%as.matrix(AntMod_mat))%*%LGain)
LGain_FixMod = 1/(t(LGain)%*%solve(t(FixMod_mat)%*%as.matrix(FixMod_mat))%*%LGain)
SGain = c(0,
0,1,-1,0,0, # small gain versus neutral
0,0,0,0,0,0,0,0,0,0
)
SGain_CueMod = 1/(t(SGain)%*%solve(t(CueMod_mat)%*%as.matrix(CueMod_mat))%*%SGain)
SGain_AntMod = 1/(t(SGain)%*%solve(t(AntMod_mat)%*%as.matrix(AntMod_mat))%*%SGain)
SGain_FixMod = 1/(t(SGain)%*%solve(t(FixMod_mat)%*%as.matrix(FixMod_mat))%*%SGain)
LGain_bl = c(0,
1,0,0,0,0, # large gain versus baseline
0,0,0,0,0,0,0,0,0,0
)
LGain_bl_CueMod = 1/(t(LGain_bl)%*%solve(t(CueMod_mat)%*%as.matrix(CueMod_mat))%*%LGain_bl)
LGain_bl_AntMod = 1/(t(LGain_bl)%*%solve(t(AntMod_mat)%*%as.matrix(AntMod_mat))%*%LGain_bl)
LGain_bl_FixMod = 1/(t(LGain_bl)%*%solve(t(FixMod_mat)%*%as.matrix(FixMod_mat))%*%LGain_bl)
SGain_bl = c(0,
0,1,0,0,0, # small gain versus baseline
0,0,0,0,0,0,0,0,0,0
)
SGain_bl_CueMod = 1/(t(SGain_bl)%*%solve(t(CueMod_mat)%*%as.matrix(CueMod_mat))%*%SGain_bl)
SGain_bl_AntMod = 1/(t(SGain_bl)%*%solve(t(AntMod_mat)%*%as.matrix(AntMod_mat))%*%SGain_bl)
SGain_bl_FixMod = 1/(t(SGain_bl)%*%solve(t(FixMod_mat)%*%as.matrix(FixMod_mat))%*%SGain_bl)
# VIF w/o intercept
#diag(solve(cor(CueMod_mat[,2:16])))
#diag(solve(cor(FixMod_mat[,2:16])))
#diag(solve(cor(AntMod_mat[,2:16])))
sub_details = c(sub, run,
LGain_CueMod, LGain_AntMod, LGain_FixMod,
SGain_CueMod, SGain_AntMod, SGain_FixMod,
# implicit baseline
LGain_bl_CueMod, LGain_bl_AntMod, LGain_bl_FixMod,
SGain_bl_CueMod, SGain_bl_AntMod, SGain_bl_FixMod)
df = rbind(df,sub_details)
}Pltting the BOLD timeseries based on each cue type from the provided behavioral data. By visualizing the different models, can observe how duration and amplitude shift between them. As a reminder, the models include:
design_CueMod: Cue Onset + Cue Durationdesign_AntMod: Cue Onset + (Cue Duration + Fixation
Duration)design_FixMod: Fixation Onset + Fixation Durationpal <- c("#000000","#004949","#009292","#ff6db6","#ffb6db",
"#490092","#006ddb","#b66dff","#6db6ff","#b6dbff",
"#920000","#924900","#db6d00","#24ff24","#ffff6d")
cuemod_plt = design_CueMod %>%
gather("Cue","Value", BW:SLmss) %>%
ggplot(aes(x = TR, y = Value, colour = Cue)) +
geom_line() +
facet_wrap(~Cue) +
ggtitle("A") +
#scale_color_brewer(palette = 'Dark2') +
scale_color_manual(values = pal) +
theme_minimal()+
theme(text = element_text(family = "Times New Roman"),
axis.text = element_text(size = 10),
axis.title = element_text(size = 10),
legend.text = element_text(size = 10),
legend.title = element_text(size = 10),
plot.title = element_text(size = 16))
antmod_plt = design_AntMod %>%
gather("Cue","Value", BW:SLmss) %>%
ggplot(aes(x = TR, y = Value, colour = Cue)) +
geom_line() +
facet_wrap(~Cue) +
ggtitle("B") +
#scale_color_brewer(palette = 'Dark2') +
scale_color_manual(values = pal) +
theme_minimal()+
theme(text = element_text(family = "Times New Roman"),
axis.text = element_text(size = 10),
axis.title = element_text(size = 10),
legend.text = element_text(size = 10),
legend.title = element_text(size = 10),
plot.title = element_text(size = 16))
fixmod_plt = design_FixMod %>%
gather("Cue","Value", BW:SLmss) %>%
ggplot(aes(x = TR, y = Value, colour = Cue)) +
geom_line() +
facet_wrap(~Cue) +
ggtitle("C") +
#scale_color_brewer(palette = 'Dark2') +
scale_color_manual(values = pal) +
theme_minimal()+
theme(text = element_text(family = "Times New Roman"),
axis.text = element_text(size = 10),
axis.title = element_text(size = 10),
legend.text = element_text(size = 10),
legend.title = element_text(size = 10),
plot.title = element_text(size = 16))
cuemod_pltantmod_pltfixmod_pltThe below is harder to visualize because the fast event design is shorter than the BOLD signal. Hence, events have quite a bit of overlap.
pal <- c("#000000","#004949","#009292","#ff6db6","#ffb6db",
"#490092","#006ddb","#b66dff","#6db6ff","#b6dbff",
"#920000","#924900","#db6d00","#24ff24","#ffff6d")
design_CueMod %>%
gather("Cue","Value", BW:SLmss) %>%
ggplot(aes(x = TR, y = Value, colour = Cue)) +
geom_line() +
ggtitle("Cue Onset + Cue Dur") +
scale_color_manual(values = pal) +
theme_minimal() +
theme(text = element_text(family = "Times New Roman"),
axis.text = element_text(size = 10),
axis.title = element_text(size = 10),
legend.text = element_text(size = 10),
legend.title = element_text(size = 10),
plot.title = element_text(size = 16))design_AntMod %>%
gather("Cue","Value", BW:SLmss) %>%
ggplot(aes(x = TR, y = Value, colour = Cue)) +
geom_line() +
ggtitle("Cue Onset + (Cue Dur + Fix Dur)") +
scale_color_manual(values = pal) +
theme_minimal() +
theme(text = element_text(family = "Times New Roman"),
axis.text = element_text(size = 10),
axis.title = element_text(size = 10),
legend.text = element_text(size = 10),
legend.title = element_text(size = 10),
plot.title = element_text(size = 16)) design_FixMod %>%
gather("Cue","Value", BW:SLmss) %>%
ggplot(aes(x = TR, y = Value, colour = Cue)) +
geom_line() +
ggtitle("Fix Onset + Fix Dur") +
scale_color_manual(values = pal) +
theme_minimal() +
theme(text = element_text(family = "Times New Roman"),
axis.text = element_text(size = 10),
axis.title = element_text(size = 10),
legend.text = element_text(size = 10),
legend.title = element_text(size = 10),
plot.title = element_text(size = 16))Plotting each effieincy model by contrast to show the estimated design efficiency for each datapoint.
# create dataframe of details w/ long format
eff_data = as.data.frame(df)
colnames(eff_data) <- c("Sub","Run",
"LarGain_CueMod", "LarGain_AntMod", "LarGain_FixMod",
"SmGain_CueMod", "SmGain_AntMod", "SmGain_FixMod",
"LarGainvBase_CueMod", "LarGainvBase_AntMod", "LarGainvBase_FixMod",
"SmGainvBase_CueMod", "SmGainvBase_AntMod", "SmGainvBase_FixMod")
eff_data_long = eff_data %>%
gather("Model","Eff",LarGain_CueMod:SmGainvBase_FixMod) %>%
separate(col = Model, into = c("Con","Model"), sep = '_')
eff_data_long$Eff <- as.integer(eff_data_long$Eff)
# plotting
# plot by run and contrast
plot_runcon = eff_data_long %>% ggplot(aes(x = Con, y = Eff, colour = Con, fill = Con)) +
# setting up raincloud
geom_rain(alpha = .3, rain.side = 'l',
boxplot.args.pos = list(
width = 0.05, position = position_nudge(x = 0.15)),
violin.args.pos = list(
side = "l",
width = 0.7, position = position_nudge(x = 0.3))) +
# significance comparison
ggsignif::geom_signif(
comparisons = list(c("LarGain", "LarGainvBase"),
c("LarGainvBase", "SmGain"),
c("SmGain","LarGain"),
c("SmGainvBase","LarGain")
),
map_signif_level = TRUE,
y_position = c(34,31,28,25,22),
margin_top = 0.1) +
# plot elements
theme_classic() +
scale_fill_brewer(palette = 'Dark2') +
scale_color_brewer(palette = 'Dark2') +
guides(fill = 'none', color = 'none') +
labs(title = "Estimated Efficiency by Run & Contrast",
caption = "LGain: Large Gain > Neut; SGain: Small Gain > Neut;
LGain v BL: Large Gain > Implicit Baseline; SGain v BL: Small Gain > Implicit Baseline")+
ylab("Estimated Effieincy")+
scale_x_discrete(name = "Model Type",
labels = c("LGain","SGain","LGain v BL","SGain v BL"))+
facet_wrap(~Run + Model)+
theme(text = element_text(family = "Times New Roman"),
axis.text = element_text(size = 10),
axis.title = element_text(size = 10),
legend.text = element_text(size = 10),
legend.title = element_text(size = 10),
plot.title = element_text(size = 16))
plot_runcon